/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://github.com/shenyczz/CSharpKit
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/

using System;
using CSharpKit.Extensions;

namespace CSharpKit.Numerics.Utils
{
    /// <summary>
    /// 拟合
    /// </summary>
    public static class KFitting
    {
        /// <summary>
        /// 名称: 最小二乘曲线拟合 <br/>
        /// 公式: y = a0 + a1*(x-xavg) + a2*(x-xavg)^2 + ...
        /// </summary>
        /// <param name="x">一维数组，长度为n，存放给定n个数据点的x值</param>
        /// <param name="y">一维数组，长度为n，存放给定n个数据点的y值</param>
        /// <param name="m">是多项式的阶数。m &lt;= n 且 m &lt;= 20；返回时,存放拟合多项式的m个系数(m0,m1,...,m-1)</param>
        /// <param name="dt">
        /// 一维数组，长度为4; 返回时,<br/>
        /// dt[0]保存拟合多项式与数据点的误差平方和<br/>
        /// dt[1]保存拟合多项式与数据点的误差绝对值之和<br/>
        /// dt[2]保存拟合多项式与数据点的误差绝对值最大值<br/>
        /// dt[3]保存x的平均值<br/>
        /// </param>
        /// <returns>a[m]</returns>
        public static double[] LeastSquares(double[] x, double[] y, int m, out double[] dt)
        {
            dt = new double[4];

            int n = x.Length;
            if (n != y.Length) throw new Exception("");

            // 要求m<=n且m<=20。若m>n或m>20，
            // 本函数自动按m=min{n,20}处理。
            if (m > n) m = n;
            if (m > 20) m = 20;
            var a = new double[m];

            double z = x.Mean();
            dt[3] = z;

            var b = new double[20];
            b[0] = 1.0;

            double d1 = n;
            double p = 0, c = 0;
            for (int i = 0; i < n; i++)
            {
                p = p + (x[i] - z);
                c = c + y[i];
            }
            c /= d1;
            p /= d1;

            a[0] = c * b[0];

            double g = 0.0, q = 0.0, d2 = 0.0;
            var s = new double[20];
            var t = new double[20];

            if (m > 1)
            {
                t[1] = 1.0;
                t[0] = -p;
                d2 = 0.0;
                c = 0.0;
                g = 0.0;
                for (int i = 0; i < n; i++)
                {
                    q = x[i] - z - p; d2 = d2 + q * q;
                    c = c + y[i] * q;
                    g = g + (x[i] - z) * q * q;
                }
                c = c / d2;
                p = g / d2;
                q = d2 / d1;
                d1 = d2;
                a[1] = c * t[1];
                a[0] = c * t[0] + a[0];
            }

            for (int j = 2; j < m; j++)
            {
                s[j] = t[j - 1];
                s[j - 1] = -p * t[j - 1] + t[j - 2];
                if (j >= 3)
                    for (int k = j - 2; k >= 1; k--) s[k] = -p * t[k] + t[k - 1] - q * b[k];
                s[0] = -p * t[0] - q * b[0];
                d2 = 0.0;
                c = 0.0;
                g = 0.0;
                for (int i = 0; i < n; i++)
                {
                    q = s[j];
                    for (int k = j - 1; k >= 0; k--) q = q * (x[i] - z) + s[k];
                    d2 = d2 + q * q; c = c + y[i] * q;
                    g = g + (x[i] - z) * q * q;
                }
                c = c / d2;
                p = g / d2;
                q = d2 / d1;
                d1 = d2;
                a[j] = c * s[j]; t[j] = s[j];
                for (int k = j - 1; k >= 0; k--)
                {
                    a[k] = c * s[k] + a[k];
                    b[k] = t[k]; t[k] = s[k];
                }
            }

            dt[0] = 0.0;
            dt[1] = 0.0;
            dt[2] = 0.0;

            for (int i = 0; i < n; i++)
            {
                q = a[m - 1];
                for (int k = m - 2; k >= 0; k--)
                {
                    q = a[k] + q * (x[i] - z);
                }
                p = q - y[i];
                if (System.Math.Abs(p) > dt[2])
                {
                    dt[2] = System.Math.Abs(p);
                }
                dt[0] = dt[0] + p * p;
                dt[1] = dt[1] + System.Math.Abs(p);
            }

            return a;
        }


        /// <summary>
        /// 切比雪夫曲线拟合
        /// y = a0 + a1*x + a2*x^2 + ...
        /// </summary>
        /// <param name="x">一维数组，长度为n，存放给定n个数据点的x值,其中x[0] &lt; x[1] &lt; ... &lt; x[n-1]</param>
        /// <param name="y">一维数组，长度为n，存放给定n个数据点的y值</param>
        /// <param name="m">是多项式的阶数。m &lt;= n 且 m &lt;= 20；返回时,存放拟合多项式的m个系数(m0,m1,...,m-1)</param>
        /// <param name="sdmax">拟合多项式偏差的最大值</param>
        /// <returns></returns>
        public static double[] Chebyshev(double[] x, double[] y, int m, out double sdmax)
        {
            int n = x.Length;
            if (n != y.Length) throw new Exception("");

            // 要求m<=n且m<=20。若m>n或m>20，
            // 本函数自动按m=min{n,20}处理。
            if (m > n) m = n;
            if (m > 20) m = 20;
            var a = new double[m];

            var ix = new int[21];
            var h = new double[21];
            var aa = new double[m + 1];

            if (m >= n) m = n - 1;
            if (m >= 20) m = 19;

            int m1 = m + 1;
            int l = (n - 1) / m;
            int j0 = l;

            ix[0] = 0;
            ix[m] = n - 1;
            for (int i = 1; i < m; i++)
            {
                ix[i] = j0;
                j0 += l;
            }

            int ii, im;
            double am, y1, y2, h1, h2, ha = 0, d, hm;

            // int itemp = 0;
            // int ltemp = 1;

            while (true)
            {
                double hh = 1.0;
                for (int i = 0; i <= m; i++)
                {
                    aa[i] = y[ix[i]];
                    h[i] = -hh;
                    hh = -hh;
                }

                for (int j = 1; j <= m; j++)
                {
                    ii = m1;
                    y2 = aa[ii - 1];
                    h2 = h[ii - 1];
                    for (int i = j; i <= m; i++)
                    {
                        d = x[ix[ii - 1]] - x[ix[m1 - i - 1]];
                        y1 = aa[m - i + j - 1];
                        h1 = h[m - i + j - 1];
                        aa[ii - 1] = (y2 - y1) / d;
                        h[ii - 1] = (h2 - h1) / d;
                        ii = m - i + j;
                        y2 = y1;
                        h2 = h1;
                    }
                }

                hh = -aa[m] / h[m];
                for (int i = 0; i <= m; i++)
                {
                    aa[i] = aa[i] + h[i] * hh;
                }

                for (int j = 1; j < m; j++)
                {
                    ii = m - j;
                    d = x[ix[ii - 1]];
                    y2 = aa[ii - 1];
                    for (int k = m1 - j; k <= m; k++)
                    {
                        y1 = aa[k - 1];
                        aa[ii - 1] = y2 - d * y1;
                        y2 = y1;
                        ii = k;
                    }
                }

                hm = System.Math.Abs(hh);
                if (hm < ha || System.Math.Abs(hm - ha) < 1.0e-12) // hm==ha
                {
                    aa[m] = -hm;
                    goto END;
                }

                aa[m] = hm;
                ha = hm;
                im = ix[0];
                h1 = hh;
                j0 = 0;
                for (int i = 0; i < n; i++)
                {
                    if (i == ix[j0])
                    {
                        if (j0 < m) j0 = j0 + 1;
                    }
                    else
                    {
                        h2 = aa[m - 1];
                        for (int k = m - 2; k >= 0; k--)
                        {
                            h2 = h2 * x[i] + aa[k];
                        }
                        h2 = h2 - y[i];
                        if (System.Math.Abs(h2) > hm)
                        {
                            hm = System.Math.Abs(h2);
                            h1 = h2;
                            im = i;
                        }
                    }
                }

                if (im == ix[0]) goto END;

                int itemp = 0;
                int ltemp = 1;
                while (ltemp == 1)
                {
                    ltemp = 0;
                    if (im >= ix[itemp])
                    {
                        itemp = itemp + 1;
                        if (itemp <= m) ltemp = 1;
                    }
                }
                if (itemp > m) itemp = m;
                if (itemp == (itemp / 2) * 2) h2 = -hh;
                else h2 = hh;
                if (h1 * h2 >= 0.0) ix[itemp] = im;
                else
                {
                    if (im < ix[0])
                    {
                        for (int j = m - 1; j >= 0; j--) ix[j + 1] = ix[j];
                        ix[0] = im;
                    }
                    else
                    {
                        if (im > ix[m])
                        {
                            for (int j = 1; j <= m; j++) ix[j - 1] = ix[j];
                            ix[m] = im;
                        }
                        else
                        {
                            ix[itemp - 1] = im;
                        }
                    }
                }

            }


        END:
            //---
            am = aa[m];
            for (int i = 0; i < m; i++)
            {
                a[i] = aa[i];
            }

            sdmax = am;
            return a;
        }

        /// <summary>
        /// 曲线拟合 (用五点三次平滑公式对等距点上的观测数据进行平滑)
        /// y = a + bx + cx^2 + dx^3
        /// </summary>
        /// <param name="y">等距点上的观测数据数组指针,返回时为等距点上的观测数据的平滑结果</param>
        public static void Curvefitting(double[] y)
        {
            int n = y.Length;
            if (n < 5) throw new Exception("");
            // 
            var yy = new double[n];

            yy[0] = 69.0 * y[0] + 4.0 * y[1] - 6.0 * y[2] + 4.0 * y[3] - y[4];
            yy[0] = yy[0] / 70.0;
            yy[1] = 2.0 * y[0] + 27.0 * y[1] + 12.0 * y[2] - 8.0 * y[3];
            yy[1] = (yy[1] + 2.0 * y[4]) / 35.0;

            for (int i = 2; i < n - 2; i++)
            {
                yy[i] = -3.0 * y[i - 2] + 12.0 * y[i - 1] + 17.0 * y[i];
                yy[i] = (yy[i] + 12.0 * y[i + 1] - 3.0 * y[i + 2]) / 35.0;
            }

            yy[n - 2] = 2.0 * y[n - 5] - 8.0 * y[n - 4] + 12.0 * y[n - 3];
            yy[n - 2] = (yy[n - 2] + 27.0 * y[n - 2] + 2.0 * y[n - 1]) / 35.0;
            yy[n - 1] = -y[n - 5] + 4.0 * y[n - 4] - 6.0 * y[n - 3];
            yy[n - 1] = (yy[n - 1] + 4.0 * y[n - 2] + 69.0 * y[n - 1]) / 70.0;

            Buffer.BlockCopy(yy, 0, y, 0, n);

            return;
        }


        // {{END}}
    }









}
